function [fX,AA] = bpass(X,pl,pu,root,drift,ifilt,nfix,thet)
%
%	MATLAB COMMAND FOR DEFAULT FILTER:  fX = bpass(X,pl,pu)
%
%	Required Inputs:
%  	X		- matrix of data (T x # variables) with variables in columns
%		pl		- minimum period of oscillation of desired component 
%		pu		- maximum period of oscillation of desired component (2<=pl<pu<infinity).
%
%  Output:
%		fX 	- matrix (T x nvars) containing filtered data 
%
%	Examples:  
%   Quarterly data: pl=6, pu=32 returns component with periods between 1.5 and 8 yrs.
%	 Monthly data  : pl=2, pu=24 returns component with all periods less than 2 yrs.
%
%  Note:  When feasible, we recommend dropping 2 years of data from the beginning 
%		and end of the filtered data series.  These data are relatively poorly estimated.
%
%	Optional Inputs:  root,drift,ifilt,nfix,thet.  THESE INPUTS ARE NOT REQUIRED.
%
%	Additional Output:  AA - T x T matrix which premultiplies a T x Nvars data matrix 
%				to compute the band-pass filtered data.  (Note: data should have any drift
%				term or time trend removed before applying AA.)
%
%=======================================================================================
%		The next discussion is only for users who wish to use optimal band pass filters 
%		other than the default filter recommended in Christiano and Fitzgerald (1999). 
%=======================================================================================
%
%  Version Date: 6/10/99 (Please notify Terry Fitzgerald at
%			tj.fitzgerald@clev.frb.org or (216)579-2293 if you encounter bugs).  
%
%  This is a matlab function that filters time series data using  
%  approximations to the band pass filter as discussed in the paper
%  "The Band Pass Filter" by Lawrence J. Christiano and Terry J. Fitzgerald (1999).

%	Several band-pass approximation strategies can be selected in this 
%  program.  To explain how to select these alternatives, we rely upon the notation 
%  developed in the paper.  The default setting used above returns the filtered data fX 
%  associated with the unrestricted optimal filter assuming a Random Walk with drift
%  time series process.  
%
%	Alternative filters require selecting the type of filter (asymetric, symmetric, or 
%  fixed length symmetric) and the underlying time series assumption (random walk, IID, 
%  or a user provided MA process with or without a unit root).  The fixed length 
%  symmetric filter advocated by Baxter and King(1999), and the trigonometric regression
%	approach, can also be chosen.  
%
%  Optional Inputs:
%		root		= 0	no unit root in time series						(default: root = 1);
%					= 1	unit root in time series
%		drift		= 0	no drift or time trend in time series			(default: drift = 1);
%					= 1	time series is Random Walk with possibly a nonzero drift
%					 		or time series is stationary about a linear time trend.
%		ifilt 	= 0	Asymmetric Filter 									(default: ifilt = 0);
%					= 1	Symmetric Filter
%					= 2	Fixed Length Symmetric Filter
%					= 3 	Baxter-King Fixed Length Symmetric Filter 
%					= 4	Trigonometric Regression Filter
%		nfix	= sets fixed length (used w/ ifilt=2 or 3)				(default: nfix=-1).
%		thet	= MA coefficients for time series model 					(default: thet=1)
%				if root=1 :	x(t) = mu + x(t-1) + thet(1)e(t) + thet(2)e(t-1) + ...
%				if root=0 :	x(t) = thet(1)e(t) + thet(2)e(t-1) + ...
%						Examples:
%									root=1, thet=1 	Random Walk filter
%									root=0, thet=1		IID filter
%
%  The following function calls are permitted:
%		fX = bpass(X,pl,pu)										: optimal Random Walk filter
%		fX = bpass(X,pl,pu,root)								: choose unit root option
%		fX = bpass(X,pl,pu,root,drift)						: choose drift/trend option
%		fX = bpass(X,pl,pu,root,drift,ifilt)				: choose type of filter
%		fX = bpass(X,pl,pu,root,drift,ifilt,nfix)			: choose fixed filter length
%		fX = bpass(X,pl,pu,root,drift,ifilt,nfix,thet)	: choose MA time series model
%
%  Notes:  
%	1) Users wishing to use a fixed length symmetric filter (ifilt = 2 or 3) must
%     also supply the lead/lag length of the filter (nfix).
%	2) The filtered data matrix (fX)	associated with some filters include zeros at the 
%		beginning and end of the data set.  The filter is not available for these points.
%		You should truncate these points from your analysis.  Examples include
%		fixed length filters and filters for user supplied time series with length(thet)>1.  
%		The default filter returns all nonzero data.
%
if nargin < 4
   root = 1;
   drift = 1;
   ifilt = 0;
   nfix = -1;
   thet = 1;
elseif nargin < 5
   drift = 1;
   ifilt = 0;
   nfix = -1;
   thet = 1;
elseif nargin < 6
   ifilt = 0;
   nfix = -1;
   thet = 1;
elseif nargin < 7
   nfix = -1;
   thet = 1;
elseif nargin < 8
   thet = 1;
end

nq=length(thet)-1;
b=2*pi/pl;
a=2*pi/pu;

[T,nvars] = size(X);

if T < 5
   warning(' (bpass): T < 5')
end
if T < 2*nq+1
   error(' (bpass): T must be at least 2*q+1');
end
if pu <= pl
   error(' (bpass): pu must be larger than pl')
end
if pl < 2
   warning(' (bpass): pl less than 2 , reset to 2')
   pl = 2;
end
if root ~= 0 & root ~= 1
   error(' (bpass): root must be 0 or 1')
end
if drift<0 | drift > 1
   error(' (bpass): drift must be 0 or 1')
end
if ifilt<0 | ifilt > 4
   error(' (bpass): ifilt must be 0, 1, 2, 3, or 3')
end
if (ifilt == 2 | ifilt == 3) & nfix < 1
   error(' (bpass): fixed lag length must be >= 1')
end
if ifilt == 2 & nfix < nq
   error(' (bpass): fixed lag length must be >= q')
end
if (ifilt == 2 | ifilt ==3) & nfix >= T/2
      error(' (bpass): fixed lag length must be < T/2')
end
if (ifilt == 4 & T-2*floor(T/2) ~= 0)
   error(' (bpass): trigonometric regressions only available for even T')
end

[m1,m2]=size(thet);
if m1 > m2 
   th=thet;
else
   th=thet';
end
%   compute g(thet)
%   [g(1) g(2) .... g(2*nq+1)] correspond to [c(q),c(q-1),...,c(1),
%                                        c(0),c(1),...,c(q-1),c(q)]
%   cc = [c(0),c(1),...,c(q)]
thp=flipud(th);
g=conv(th,thp);
cc = g(nq+1:2*nq+1);
%   compute "ideal" Bs
j = 1:2*T;
B = [(b-a)/pi (sin(j*b)-sin(j*a))./(j*pi)]';
%    compute R using closed form integral solution
R = zeros(T,1);
if nq > 0
   R0 = B(1)*cc(1) + 2*B(2:nq+1)'*cc(2:nq+1);
   R(1) = pi*R0;
   for i = 2:T
      dj = Bge(i-2,nq,B,cc);
      R(i) = R(i-1) - dj;
   end
else
   R0 = B(1)*cc(1);
   R(1) = pi*R0;
   for j = 2:T
      dj = 2*pi*B(j-1)*cc(1);
      R(j) = R(j-1) - dj;
   end
end
%
%======================================================================
%
AA = zeros(T,T);
%
%----------------------------------
%  asymmetric filter
%----------------------------------
%
if ifilt == 0
   if nq==0
      for i = 1:T
         AA(i,i:T) = B(1:T-i+1)';
         if root == 1;
            AA(i,T) = R(T+1-i)/(2*pi);
         end
      end
      AA(1,1) = AA(T,T);
      %  Use symmetry to construct bottom 'half' of AA
      AA = AA + flipud(fliplr(triu(AA,1)));
   else
      % CONSTRUCT THE A MATRIX size T x T
      A = Abuild(T,nq,g,root);
      Ainv = inv(A);
      % CONSTRUCT THE d MATRIX size T x 1 
      for np = 0:ceil(T/2-1)
         d = zeros(T,1);
         ii = 0;
        	for jj = np-root:-1:np+1+root-T
           	ii = ii+1;
	         d(ii) = Bge(jj,nq,B,cc);
   	   end
         if root == 1
            d(T-1) = R(T-np);
         end
         %  COMPUTE Bhat = inv(A)*d
         Bhat = Ainv*d;
         AA(np+1,:) = Bhat';
      end
      %  Use symmetry to construct bottom 'half' of AA
      AA(ceil(T/2)+1:T,:) = flipud(fliplr(AA(1:floor(T/2),:)));
   end
end   
%
%-------------------------------------
%  symmetric filter
%-------------------------------------
%
if ifilt == 1
   if nq==0
      for i = 2:ceil(T/2)
         np = i-1;
         AA(i,i:i+np) = B(1:1+np)';
         if root == 1;
            AA(i,i+np) = R(np+1)/(2*pi);
         end
         AA(i,i-1:-1:i-np) = AA(i,i+1:i+np);
      end
      %  Use symmetry to construct bottom 'half' of AA
      AA(ceil(T/2)+1:T,:) = flipud(fliplr(AA(1:floor(T/2),:)));
   else
      for np = nq:ceil(T/2-1)
         nf = np;
         nn = 2*np+1;
      % CONSTRUCT THE A MATRIX size nn x nn
       	A = Abuild(nn,nq,g,root);
         Ainv = inv(A);
      % CONSTRUCT THE d MATRIX size nn x 1
         d = zeros(nn,1);
         ii = 0;
         for jj = np-root:-1:-nf+root
  	         ii = ii+1;
     	      d(ii) = Bge(jj,nq,B,cc);
         end
         if root == 1
            d(nn-1) = R(nf+1);
			end            
      %  COMPUTE Bhat = inv(A)*d
         Bhat = Ainv*d;
         AA(np+1,1:2*np+1) = Bhat';
      end
      %  Use symmetry to construct bottom 'half' of AA
      AA(ceil(T/2)+1:T,:) = flipud(fliplr(AA(1:floor(T/2),:)));
   end
end
%
%----------------------------------
%  fixed length symmetric filter
%----------------------------------
%
if ifilt == 2
   if nq==0
      bb = zeros(2*nfix+1,1);
      bb(nfix+1:2*nfix+1) = B(1:nfix+1);
      bb(nfix:-1:1) = B(2:nfix+1);
      if root == 1
         bb(2*nfix+1) = R(nfix+1)/(2*pi);
         bb(1) = R(nfix+1)/(2*pi);
      end
      for i = nfix+1:T-nfix
         AA(i,i-nfix:i+nfix) = bb';
      end
   else
      nn = 2*nfix+1;
      % CONSTRUCT THE A MATRIX size nn x nn
     	A = Abuild(nn,nq,g,root);
      Ainv = inv(A);
      % CONSTRUCT THE d MATRIX size nn x 1
      d = zeros(nn,1);
      ii = 0;
      for jj = nfix-root:-1:-nfix+root
 	      ii = ii+1;
	  	   d(ii) = Bge(jj,nq,B,cc);
      end
      if root == 1
     	   d(nn-1) = R(nn-nfix);
      end
      %  COMPUTE Bhat = inv(A)*d
      Bhat = Ainv*d;
      for ii = nfix+1:T-nfix
         AA(ii,ii-nfix:ii+nfix) = Bhat';
      end
   end
end
%
%----------------------------------
%  Baxter-King filter
%----------------------------------
%
if ifilt == 3
   bb = zeros(2*nfix+1,1);
   bb(nfix+1:2*nfix+1) = B(1:nfix+1);
   bb(nfix:-1:1) = B(2:nfix+1);
   bb = bb - sum(bb)./(2*nfix+1);
   for i = nfix+1:T-nfix
      AA(i,i-nfix:i+nfix) = bb';
   end
end
%
%----------------------------------
%  Trigonometric Regression filter
%----------------------------------
%
if ifilt == 4
	jj = 1:T./2;
	% find frequencies in desired band omitting T/2;
	jj = find(T./pu<=jj & jj<=T./pl & jj<T./2);
	if isempty(jj)
	   error(' (bpass): frequency band is empty in trigonometric regression');
	end
	om = 2*pi*jj./T;
	if pl > 2
  		for t = 1:T
   		for k = T:-1:1
      		l = t-k;
      		tmp = sum(cos(om*l));
      		AA(t,k) = tmp;
   		end
  		end
	else
  		for t = 1:T
   		for k = T:-1:1
      		l = t-k;
      		tmp = sum(cos(om*l));
      		tmp2 = ( cos(pi*(t-l)).*cos(pi*t) )./ 2;
      		AA(t,k) = tmp + tmp2;
   		end
  		end
	end
	AA = AA*2/T;
end
%
%======================================================================
%
%  check that sum of all filters equal 0 if assuming unit root
%
if root == 1
   tst = max(abs(sum(AA')));
   if tst > 1.e-09 & root ~= 0
      warning(' (bpass): Bhat does not sum to 0 ')
      tst   
   end
end
%
%======================================================================
%
%	compute filtered time series using selected filter matrix AA
%
if drift > 0;
   X = undrift(X);
end
fX = AA*X;
%
%======================================================================
%				Functions
%======================================================================
%
function y = Bge(jj,nq,B,cc)
%
%  closed form solution for integral of B(z)g(z)(1/z)^j  (eqn 16)
%     nq > 0, jj >= 0
%     if nq = 0, y = 2*pi*B(absj+1)*cc(1);
%
absj =abs(jj);
if absj >= nq
   dj = B(absj+1)*cc(1) + B(absj+2:absj+nq+1)'*cc(2:nq+1);
   dj = dj + flipud(B(absj-nq+1:absj))'*cc(2:nq+1);
elseif absj >= 1
   dj = B(absj+1)*cc(1) + B(absj+2:absj+nq+1)'*cc(2:nq+1);
   dj = dj + flipud(B(1:absj))'*cc(2:absj+1);
   dj = dj + B(2:nq-absj+1)'*cc(absj+2:nq+1);
else
   dj = B(absj+1)*cc(1) + 2*B(2:nq+1)'*cc(2:nq+1);
end
y = 2*pi*dj;
%
%-----------------------------------------------------------------------
%
function A = Abuild(nn,nq,g,root)
%
%  builds the nn x nn A matrix in A.12
%   if root == 1 (unit root)
%      Abig is used to construct all but the last 2 rows of the A matrix
%   elseif root == 0 (no unit root)
%      Abig is used to construct the entire A matrix
%
if root == 1
	Abig=zeros(nn-2,nn+2*(nq-1));
	for j = 1:nn-2      
	   Abig(j,[j:j+2*nq]) = g'; 
	end
	A = Abig(:,[nq:nn+nq-1]);
	%   construct A(-f)
	Q = -ones(nn-1,nn);
	Q = tril(Q);
	F = zeros(1,nn-1);
	F(nn-1-nq:nn-1) = g(1:nq+1);
	A(nn-1,:) = F*Q;
	%    construct last row of A
   A(nn,:) = ones(1,nn);
else
	Abig=zeros(nn,nn+2*(nq-1));
	for j = 1:nn      
   	Abig(j,[j:j+2*nq]) = g'; 
	end
	A = Abig(:,[nq+1:nn+nq]);
end   
%    multiply A by 2*pi
A = 2*pi*A;
%
%-----------------------------------------------------------------------
%
function xun = undrift(x)
%
%  This function removes the drift or a linear time trend from a time series using the formula
%			drift = (x(T) - x(1)) / (T-1).
%
%  Input:  x - data matrix x where columns represent different variables, x is (T x # variables). 
%  Output: xun - data matrix same size as x with a different drift/trend removed from each variable.
%
[T,nvars] = size(x);
xun = zeros(T,nvars);
dd = [0:T-1]';
for ivar = 1:nvars
   drift = (x(T,ivar)-x(1,ivar)) / (T-1);
   xun(:,ivar) = x(:,ivar) - dd*drift;
end



